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Abstract 

Scaling properties of patterns formed by large contact forces are studied as a function of the 
applied shear stress, in two-dimensional static packings generated from the force network ensemble. 
An anisotropic finite-size-scaling analysis shows that the applied shear does not affect the universal 
scaling properties of these patterns, but simply induces different length scales in the principal 
directions of the macroscopic stress tensor. The ratio of these length scales quantifies the anisotropy 
of the force networks, and is found not to depend on the details of the underlying contact network, 
in contrast with other properties such as the yield stress. 

PACS numbers: 45.70.-n, 45.70.Cc, 46.65. +g 
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Aggregates of macroscopic particles such as granular materials, foams and emulsions are 
often found in a disordered, solid-like state whose mechanic al p roperties have attracted much 



attention in recent years 
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131 ]. Under external stresses, 



these systems present a non-zero yield threshold which is solely due to the intricate network 
formed by contact forces between particles. One of the most remarked features of these 
highly disordered force networks is the tendency for large forces to align and form branching 
force chains. When these appear in response to an external stress, such as a global shear, 
their spatial inhomogeneity is most striking jf], 3, Is]- As the applied shear stress is increased, 
the large forces orient in a preferred direction, and the force network becomes more and more 
anisotropic (c.f. Fig. [Q, up to the point where the applied stress can no longer be sustained, 
and the packing yields. While this qualitative picture has been long established, due to 
a lack of appropriate analytic and numerical tools, quantitative studies of packings under 
static shear have been few 



Recently, a novel characterization has been introduced for the geometrical patterns formed 
by large forces in isotropically compressed force networks 13|])- The patterns, displayed 
by molecular dynamics (MD) simulations of granular packings, turned out to have scaling 
properties which are independent of the pressure, polydispersity and friction. An intriguing 
question is whether an analysis of scaling properties can provide insight into the organization 
of force networks in packings under shear. This is a much more difficult problem. 

Creating static packings under shear with MD is seriously hindered by local rearrange- 
ments that seem to prevent the system from reaching a clear mechanical equilibrium. To 
allow for systematic examinations of the effects of shear, it has been proposed 10, LLI[ to 
ignore the microscopic rearrangements of the grains, and study the ensembles of force net- 
works allowed on a given contact network, as function of macroscopic stresses. This purely 
statistical approach was found to account well for the properties of packings under shear 
such as the existence of a yield threshold [111 ], as well as for the scaling properties of clus- 
ters of large forces in packings without shear 13J]. It also describes remarkably well the 



distribution of force magnitudes 
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121 ] and the response to an external overload 



In this Letter, we study the geometry of patterns of large forces in two-dimensional, static 
packings under shear generated from the force network ensemble. We find that applying 
shear stress leaves unchanged the universal scaling properties of the patterns, but introduces 
two different length scales in the two directions of principal stress axes. We characterize the 
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FIG. 1: Effect of shear on force networks: (a-c) Force networks obtained from the force ensemble, 
for three different values of shear: (a) r = 0, (b) r = 0.2, (c) r = 0.4. The forces are represented 
as bonds connecting centers of grains of contact, the width of each bond being proportional to 
the magnitude of the force. The underlying contact network is the same for all values of r. (d-f) 
Clusters of forces larger than a threshold / (bold bonds), obtained for / close to the critical value. 

anisotropy of the force networks by the ratio r of these length scales. An anisotropic finite- 
size scaling analysis allows us to determine r as function of the applied shear stress r and 
the coordination number of the packing. In contrast with other properties such as the 
yield stress, r(r) turns out to be independent of the underlying contact geometry, and thus 
provides a universal characterization of shear-induced anisotropy. 

Force networks — To study systematically the influence of the shear stress on the geom- 
etry of large forces we examined force networks generated from the force network ensemble 
pertaining to packings of frictionless, non-attracting particles. This approach relies on the 
fact that in a typical packing, the number of unknown contact forces is larger than the num- 
ber of equations for mechanical balance Q. For a given, fixed, contact geometry, the force 
network ensemble consists of all force networks with repulsive contact forces in mechanical 
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balance on each grain, and consistent with the tensor cr of applied macroscopic stress. The 
considered contact networks are isotropic, and the coordinates are such that the pressure is 
isotropic, i.e. a xx = a yy , while the dimensionless shear stress r = a xy /a xx is varied. 

The ensemble is sampled numerically using two different methods : (i) the wheel-move 
algorithm 12j for ordered, hexagonal packings (z = 6); (ii) the procedure outlined in for 
disordered packings of various coordination numbers z. Both methods efficiently generate 
very large packings (10 5 and 5.10 4 particles respectively), from which smaller sub-systems 



y 



or r smaller than 



were extracted and analyzed. Note that the ensemble is non-empty on 
a maximal value T max (z), which has been identified as the yield stress 

Force chains — To stud y t he geometry of patterns formed by large forces, we follow 
the method introduced in [13 1. which we here briefly recall and then extend to packings 
under shear. The force network in a packing of frictionless particles can be represented as 
a set of bonds connecting grains in contact, where each bond carries the magnitude of the 
corresponding contact force (c.f. Fig. [T] (a-c)). A natural way to isolate the force chains 
in such a network is to choose a threshold /, and consider only the subgraph formed by 
bonds carrying forces larger than / (c.f. Fig. [1] (d-f)). Varying this threshold allows to 
study patterns of large forces at different scales: for small values, most of the grains remain 
connected, but as the threshold is increased, the extracted subgraph breaks into disconnected 
clusters. The extent of each chain of forces larger than the threshold can be characterized 
by the size of the corresponding cluster, quantified by the number s of mutually connected 
bonds. 

In an ensemble of packings, a statistical description of the fluctuating patterns of large 
forces is given by the statistics of the clusters obtained at different thresholds. Following 



the methods of percolation theory 



151 ]. the geometry of an ensemble of packings can be 



characterized by P(s, f), the number (per bond) of clusters of size s at the threshold /. For 
packings under isotropic pressure, it is moreover natural to describe the clusters by a single 
characteristic length, the cluster correlation length £(/), which corresponds to the typical 
linear size of the clusters at a given threshold (not considering clusters percolating through 
the whole system). 

Scaling analysis — Analogously to percolation and other lattice models of critical phe- 
nomena, the system is critical around the value f c of the threshold above which no infinite 
cluster is found. At this value, £(/) diverges as ]/ — f c \~ u , and P(s, f) becomes a power- law, 
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its moments (s n ) diverging as powers of £(/): (s n ) ~ £^ n . For systems of finite size, the 
correlation length is naturally bounded by the linear size of system, L, and for £ comparable 
to L, the behavior of (s 11 ) is given by a scaling function of L/£ [ly, Il7| : 

(*»(/, L^L^E^/-/^ 1 /*], (1) 

where n and z/ are universal critical exponents, and E n universal scaling functions. For 
force networks under isotropic pressure, 2 = 1.89, v = 1.65, and E 2 is independent of 
pressure, polydispersity, friction and force law [131 ] . From now on, we will only consider the 
case n = 2, and omit the indices n in the right-hand-side. Note that the scaling function E 2 
varies smoothly and displays a single maximum close to zero. 

Equation ([T|) assumes that the size of the examined domains is isotropic. In the following, 
we will need to consider rectangular domains of size Li xL 2 , in which case the scaling function 
E also depends on the aspect ratio L%jL\. Considering only the maximum of {s 2 (f, Li, L 2 )) 
with respect to /, which is equivalent to looking only at the behavior at the effective critical 
point, Eq. (TjQ) reduces to 

( a 2 (L 1 ,Z 2 )) wlflB = 2^sf^y (2) 



M, 

The scaling function E(;r) simply expresses the fact that if L\ <C L 2 , the effective correlation 
length is set by Li, so that (s 2 ) max ~ L{, and correspondingly if I 2 < (s 2 ) ma x ~ L%. 
We therefore have (with y = Li/L 2 ) 

{y^ for y < 1 
(3) 
1 for y > 1 

For L\ ~ L 2 , there is a crossover between the two asymptotic trends, and the precise behavior 
is determined by the scaling function. 

Shear-induced anisotropy — The analysis presented so far pertains to the case of isotropic 
pressure. In packings under shear, the force networks become increasingly anisotropic as the 
shear stress is increased, the large forces aligning preferentially in the direction of the max- 
imal stress axis (cf. Fig. [T]). In consequence, close to the percolation threshold the clusters 
can no longer be described by a single length scale. Instead two correlation lengths £a/ and 
£ m must be distinguished, in the direction of maximal and minimal stress axes respectively. 
In principle, at criticality these two length scales could diverge with two different exponents 
Vm and u m . A finite-size scaling analysis as in [13], but now varying independently the 



system sizes Lm and L m in the directions of principal stresses, shows that vu = v m = u, 
the exponents and v being equal to those of isotropically compressed packings. 

Although they diverge with the same exponent, the values of the correlation lengths £m 
and £ m clearly differ. This implies that, to leading order, 

6/ = - fc\- u and U = $>|/ - f c \- u - (4) 

where the length scales an d depend on the shear stress r but not on the threshold 
/. In the language of critical phenomena, this situation is called weakly anisotropic scaling, 
as opposed to the case where the two correlation-length exponents differ jl?) ]. 

To quantify the anisotropy of force networks as function of shear stress, from Eq. (j3J), it 
appears natural to use the anisotropy ratio r(r) defined as 

AO) 

'M = %■ <«) 

SAT 

Clearly, for r = 0, r = 1, and as r is increased, r decreases below one. 

The anisotropy ratio r as function of r can be determined using an anisotropic scaling 
analysis. The central observation is that, in the context of conventional critical phenomena, 
a weakly anisotropic system can be made isotropic simply by rescaling the lengths in the two 
directions of principal axes. In the present setting, this property suggests that if the actual 
lengths L M and L m are replaced by properly rescaled effective lengths L M and L m , the 
scaling properties of the sheared force networks should be identical to those of isotropically 
compressed networks. This property allows us to determine the ratio of length scales ^ 
and 

More specifically, the scaling of the maximum of (s 2 ) is described by Eq. ()2]) with an 
additional dependence on the shear stress r in the right hand side. Our hypothesis is that 
the scaling function S does not depend explicitly on r, and that the anisotropy can be 
eliminated simply by replacing Lm and L m by rescaled lengths Lm and L m given by 



b M {r)L M and L m = b m (r)L m (6) 



with \)m = Vw an d b m = l/d)?, and therefore r(r) = bM/b m - Substituting into Eq. (j2J), 
the scaling relation becomes 

(s (Lm , L m , r)) max (6 m Z/ m ) S 

TV 1 • (7) 



where the dependence on r occurs only through the scale factors b m and bu- 

The validity of the scaling relation (J7|) can be checked directly from the numerical data. 
Fig. [2] illustrates the behavior of (s 2 ) max as function of L m , for three different values of r, 
and in each case three different values of Lm- For fixed Lm and r, as L m increases, at first 
a clear power-law can be observed. For larger values of L m , a crossover occurs and (s 2 ) max 
reaches a plateau. As Lm is increased, the behavior in the scaling regime is unchanged, 
but the crossover occurs at larger L m , and the value of the plateau increases. According to 
Eq. (T5]), the value of the crossover should scale as L M , and the value of the plateau as L M . 
Rescaling both axes appropriately, the data for different Lm indeed collapse on the same 
curves as shown in Fig. [2] (b). 

For different values of the shear stress, the exponents in the scaling regimes appear to 
be identical, however the prefactors of the power laws clearly depend on r. From Eq. ([7j), 
these prefactors correspond to 6 a/ and b m , which can thus be determined respectively from 
the value reached at the plateau and the intercept of the power law. Replacing Lm and L m 
by rescaled lengths Lm and L m defined in Eq. (EJ), if the scaling function E is independent 
of r, all the data must collapse on a single curve. Note that as r is increased, the cross 
over occurs at smaller L m , and the domain of power-law scaling shrinks, so that b m can be 
extracted only for r sufficiently smaller than the yield stress r max (beyond which the force 
network ensemble is empty). 

Fig. [3] displays the results obtained for all considered values of shear stress r, on several 
packing geometries of different coordination number z (z = 6 corresponds to the ordered 
hexagonal case, while the other packings are disordered). All the data clearly collapses 
on the same curve, in agreement with the hypothesis that the scaling function does not 
depend on r. As postulated, a simple rescaling of the length scales in the two directions of 
anisotropy is thus sufficient to recover isotropic scaling in force networks under shear. The 
inset of Fig. [3] displays the anisotropy parameter r = ^ as function of r. Unexpectedly, 
r(r) appears to be independent of the coordination number and geometry of the underlying 
contact network. 

Discussion — Our results show that applying an external shear stress on a force network 
does not affect the universal scaling properties of the force chains, but only induces two 
different length scales in the directions of the two principal stress axes. While these typ- 
ical lengths and their ratio r are not a priori expected to be universal, we find that they 
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FIG. 2: (color online) Scaling of the maximum of the second moment of cluster sizes: (a) 
{s 2 )max{LM,L m ) as function of L m . Three different values of Lm are represented with circles, 
squares and triangles, while three different values of the shear stress are shown in three different 
colors. The displayed data corresponds to hexagonal packing of disks, i.e. z = 6. (b) Collapses of 
data corresponding to same r, obtained by rescaling the axes by Lm- The curves are colored in 
the order of the legend. 
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FIG. 3: (color online) The anisotropic scaling function £ (cf. Eq. ([2])): collapse of data obtained 
for different values of shear stress r and different packing geometries. Each symbol corresponds to 
data for a different coordination number z (z = 6: regular, hexagonal packing, already displayed 
in Fig. [21 all others: disordered packings). Each value of r is represented by a different color. The 
inset shows the anisotropy parameter r = ^ as function of r, for different z. The black solid line 
is a quadratic fit. 
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are identical for various contact networks we considered, which include the regular hexag- 
onal packing, and disordered networks of different coordination numbers, and in this sense 
universal. 

An important remaining question is the behavior close to the yielding point. If yielding 
is analogous to a phase transition, as suggested by the jamming picture 0, S], it could 
be expected that close to it a cross-over occurs, and the scaling properties of clusters of 
large forces change significantly. The value of the yield stress T m „ x was found to be strongly 
dependent on the coordination number z of the contact network 11( . On the other hand, we 
find that r(r) is completely independent of z up to T max (z). This observation suggests the 
absence of any diverging or vanishing length scale which would accompany the cross-over 
close to T max . Moreover, we have not observed any dramatic change in the scaling properties 
close to T max , but additional work is necessary to clarify these issues. 
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